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Abstract 

In an earlier paper, Misner (2004, Class. Quant. Grav., 21, S243) presented a novel 
algorithm for computing the spherical harmonic components of data represented on 
a cubic grid. I extend Misner’s original analysis by making detailed error estimates 
of the numerical errors accrued by the algorithm, by using symmetry arguments to 
suggest a more efficient implementation scheme, and by explaining how the algo- 
rithm can be applied efficiently on data with explicit reflection symmetries. 
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1 Introduction 


Spherical harmonic analysis of data is a common procedure in many appli- 
cations in science and engineering. In numerical calculations structured on a 
cubic grid, however, extracting spherical harmonic components can be non- 
trivial since the spherical harmonic components of a function $ are defined 
by integrals over spheres 

Qimit, r) = j> Y lm (6, <p)<P{t, r, d , <p)a ,2 Q. (1) 
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(8). For a finite volume, this weighing scheme provides an approximation that 
scales as O(h), but, for a region that is also scaling with h, the resulting integral 
scales as 0(h 2 ). This provides additional motivation for choosing A oc h since 
most applications will require at least second order accuracy. 

For higher than second order accuracy, a new scheme for computing the volume 
integrals is needed, but the rest of the analysis here holds true. Given such a 
scheme, the analysis here shows how to choose the remaining parameters to 
ensure that numerical errors scale like any desired power of the grid spacing. 


4 Choosing the parameters 


In practice, the grid spacing parameter h is usually chosen to resolve the 
sources without exceeding the physical limits of the computer. I would not 
expect, in general, that the grid spacing would be chosen based on the needs 
of this algorithm. For that reason, let me assume now that h is chosen, and 
discuss how to choose the remaining parameters N and A. In this section I will 
discuss some of the theoretical issues that should be considered when choosing 
the parameters, leading to a rule of thumb that is valid based on this analysis 
and my experience with the algorithm. 

The error analysis of Section 3 implies that for fixed A, increasing the value 
of N decreases the error term. It also implies that for fixed N, increasing 
the value of A increases the error term. This suggests taking A as small as 
possible, and N as large as possible to make the error term as small as possible. 
This must be balanced, however, against practical limitations. Certainly the 
shell thickness A needs to be large enough so that there are some grid points 
within the shell, otherwise the whole procedure is undefined. For fixed N, a 
stronger restriction requires that the Legendre polynomial P W can be resolved 
over the shell. Without this condition, there would seem to be no benefit to 
taking higher values of N. Getting higher accuracy in practice requires finding 
a proper balance between choosing A small and N large. 

In making this balance, however, one must keep in mind that the error in the 
method is partially determined by the weighting scheme (8), which is only 
second order accurate in the grid spacing. I am, in addition, going to choose 
A oc h for reasons described above. This already suggests that taking N larger 
than two is pointless, since choosing N — 2 already makes the piece of the 
error that is proportional to A scale like 0{ A 4 ) (cf. Table 1), meaning that it 
will be an error term of sub-leading order in grid spacing. But once this term is 
of sub-leading order, it is much less important how large I choose A, provided 
that I still choose it proportional to the grid spacing. I therefore adopt the 
following 
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Rule of Thumb: Choose N just large enough to ensure that the error 
term proportional to A is an error term of sub-leading order in grid spacing. 
Choose A just large enough to safely resolve Pn on the shell. 


With this rule of thumb, and the second order accurate weighting scheme (8), 
I found the choices N — 2 and A = 3/i/4 completely satisfactory for a second 
order accurate code. Note that this corresponds to Misner’s choice of A in Ref. 
[1], With N = 21 found that larger values of A are also acceptable. Numerical 
results justifying these estimates appear in Ref. [2], 


5 Symmetry issues 


There are two points of interest related to this method of spherical harmonic 
decomposition and symmetries. The first was mentioned briefly in Section 2, 
namely that symmetry causes the metric G AB to be real and sparse. The sec- 
ond deals with implementing the method for grids in which explicit symmetries 
are enforced on grid functions in order to reduce the computational load of the 
simulation. In these cases, in which data is not evolved over a whole extraction 
sphere, additional analysis is required to demonstrate that the method is well 
defined and to understand how to most efficiently implement it. The primary 
result on this second topic is that the adjoint harmonics Y A of (11) have the 
same symmetries under reflection as the original spherical harmonics Y A . 


Consider first the implications of symmetry on the metric Gab ■ The symme- 
tries of the spherical harmonics, summarized in Table 2, cause the imaginary 
part of all terms in the integral (9) to cancel in pairs of points on the sphere 
related by reflections through coordinate planes. The reason is that each of the 
four signs (+1, (— l) m , (— 1) J , and (— l) z+m ) appears twice in Table 2, once for 
a term that is complex conjugated and once for a term that is not. The matrix 
is also sparse. By similar reasoning, for certain values of l and m, the terms in 
the integral (9) can cancel in sets of four. Both of these facts can be seen at 


once through a simple calculation. The ide, 


idea is to break the integral into parts 


using the second and third columns of Table 2, and then to simplify using 
the last two columns. Considering first just the symmetries under reflection 
through the xy - plane and recalling the definition (9), 
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Planes 

0 

<P 

Sign 

Conjugate 

None 

0 
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+1 

no 
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(-l) m 
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-<t> 

+1 
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Z 

7 T - 6 

4 > 

^iy+m 

no 

xy 

0 

TT + <f> 

(— !) m 

no 

XZ 

71-0 

i 

7 T — (j) 

i 

(- 1 )' 

yes 

yz 

7 T — 0 

-<t> 


yes 

xyz 

71-9 

IT + <fi • 

(- 1 )' 

no 


Table 2 

The table shows how the arguments of spherical harmonics transform under reflec- 
tions through various Cartesian planes. The first column indicates which coordinates 
have their signs inverted, while the second and third columns give the new angular 
arguments to the spherical harmonic Y[ m . Alternatively, the fourth and fifth column 
show, respectively, the overall sign in front of and whether or not to complex conju- 
gate the given spherical harmonic with the original angular arguments. The second 
row, for example, says that Yi m (-x,y,z) = Y lm (0,ir - <j>) = (-1 ) m Yi m (0 , <f>) , where 
(0,<j>) are the angular coordinates of the point ( x,y,z ). Note that this table differs 
slightly from that in Ref. [2]. The table here is correct. 

Gab = j> Yi imi (0, <t>)Yi 2m2 (0, (j>)d 2 Q (16a) 

r2ir p7:/2 _ 

= / 

Jo Jo 

r27T rn / 2 

+ Jo Jo - 0 >^)^2m 2 (7T - 9,<fi)d 2 Q (16b) 

r2ix rn/2 

= [l+(-l)‘‘ +b ]y o l Y hm ,(6,<t,)Y hm ,(e,4,)d 2 n. (16c) 

(I have suppressed the radial functions since they play no role here.) Repeating 
the procedure for reflections through the xz- and yz-planes shows that 

Gab = 2a mi +m 2 ,h+h J Q / J Q ' Re {^! mJQ , 4>)Yi 2m2 (9 , 4>)} d 2 7l (17) 

where 

<T mi +m 2 ,h+h = 1 + + (_i)h+i 2 + (-!)”*» +^+h+/ 2 (ig) 

This proves that the matrix is real. In addition, the matrix element is zero by 
symmetry whenever 

&mi+m.2, h+l 2 ~ 0 (19) 

which is true for 56 of the 81 matrix elements that exist when considering 
a fixed value of n and all values of l and m for l < 2. Of the remaining 
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mi + m 2 

h + h 

Number 

Satisfies (19) 

even 

even 

25 

no 

even 

odd 

16 

yes 

odd 

even 

20 

yes 

odd 

odd 

20 

yes 


Table 3 

The table summarizes which entries of Gab identically vanish because of the sym- 
metries of the spherical harmonics under reflections through coordinate planes for 
all values of l and m with l < 2. (This is governed by equation (19).) Of the 81 
possible matrix elements, only 25 have non-trivial values. 


25 matrix elements, 9, of course, are the diagonal elements that go to unity 
in the continuum limit. The exact break-down of which such elements must 
be zero by symmetry is summarized in Table 3. Knowing that the matrix is 
real-symmetric and sparse allows for a more efficient implementation of the 
algorithm in general. It is also extremely useful in analyzing the algorithm in 
the context of the second topic of this section, explicit grid symmetries. 


When evolving initial data with known symmetries, it is very common to 
evolve only that part of the data that is unique. In such cases, an appropriate 
symmetry boundary condition is applied at some edges of the grid. This is, 
however, inconvenient for wave extraction since computing spherical harmonic 
components (by any method) requires integrating over the full sphere. If data 
with octant symmetry, for example, is evolved only in a single octant, it is 
neither sufficient to apply the decomposition algorithm in that one octant nor 
to multiply the result of a single octant by 8 since the symmetry may forbid 
some modes as well as repeating them. 


In principle the problem appears to be even more difficult for this particular 
decomposition method. Although the spherical harmonics have well defined 
symmetries under reflections, as summarized in Table 2, it is the adjoint har- 
monics that appear in (11). The adjoint harmonics, however, are constructed 
by contracting G AB with the usual spherical harmonics, and this appears to 
mix different values of l and m. While this mixing does occur, the the matrix 
Gab is sparse in just the right way to ensure that the adjoint harmonics have 
the same symmetries at the usual spherical harmonics. 


A particular choice of the mapping (n, /, m) h* A makes this easiest to see. 
Specifically, considering all values of l and m with l < 2, there is a basis in 
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which Gab takes block diagonal form 



with all unwritten entries identically zero by symmetry. In this expression Ex 
is a AN x AN matrix over the basis functions B\ = {Y n0 o, Y n2> _ 2 , Y n2 o, Y n22 }; 
E 2 is an TV x iV matrix over the basis functions B 2 = {Y n i 0 }; E 3 is a 2N x 2 N 
matrix over the basis functions B 3 = {Y n i _i, V^n}; and E 4 is a 2N x 2N 
matrix over the basis functions B± — {Y n 2 - 1 , Y n2 i}. In this basis the matrix 
is block diagonal, so the inverse matrix 



is also block diagonal and the different basis sectors do not mix. This last point 
is key. It implies that any particular adjoint harmonic is a linear combination 
of spherical harmonics from a single set Bk 

ynlm = £ y^, ( 2 2 ) 

Y n’l’m’€ B k 

where k is the index such that Y n i m e Bk . Because, in each set Bk, the parity 
of l and the parity of m is the same on each Y n i m € B k , and because, in light 
of Table 2, it is the parity of l and m that determines the symmetries of Y n i m 
under reflections through planes, every spherical harmonic in B k for any fixed 
k has the same symmetries under reflections as any other spherical harmonic 
in B k . This implies that the adjoint harmonics also share this symmetry under 
reflection. 


6 Discussion 

In this paper I have provided detailed error estimates of the Misner algorithm 
for computing spherical harmonic components of data represented on a cubic 
grid. This analysis allows one, in principle, to chose the parameters of the 
algorithm such that its numerical errors scale any desired power of the grid 
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spacing. The only limitation of this in practice is finding a scheme for approx- 
imating volume integrals on a shell of sufficiently high accuracy. (Misner’s 
; original choice allows for a second order accurate result.) 

In addition, analysis of the symmetry properties of the spherical harmonics 
provides two insights: First, the number of operations required to initialize 
the data structures used to compute spherical harmonic components can be 
reduced by computing only those elements of Gab that are not forbidden by 
symmetries, and, second, that the adjoint harmonics used by the algorithm 
have the same symmetries under reflections as the usual spherical harmonics. 
This second fact allows the method to be efficiently used on data with explicit 
grid symmetries when only the independent portion of that data is evolved. 
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